Mon. Not. R. Astron. Soc. 000, 000-000 (0000) Printed 1 February 2008 (MN MfeK style file vl.4) 



Probing CMB Non-Gaussianity Using Local Curvature 



Olivier Dore, Stephane Colombi and Frangois R. Bouchet 

Institut d'Astrophysique de Paris, 98bis boulevard Arago, 75014 Paris, FRANCE 
dore@iap.fr, colombi@iap.fr, bouchet@iap.fr 



1 February 2008 



(N 
O 
O 
(N 

X> 
<D 
IX, 

> 

in 

m 

(N 
O 
(N 
O 

^ 

Oh 
6 

03 



X 



ABSTRACT 

It is possible to classify pixels of a smoothed cosmic microwave background (CMB) 
fluctuation map according to their local curvature in "hill" , "lake" and "saddle" re- 
gions. In the Gaussian case, fractional areas occupied by pixels of each kind can be 
computed analytically for families of excursion sets as functions of threshold and mo- 
ments of the fluctuation power spectrum. We show how the shape of these functions 
can be used to constrain accurately the level of non-Gaussianity in the data by ap- 
plying these new statistics to an hypothetical mixed model suggested by Bouchet et 
al. (2001). According to our simple test, with only one 12.5 x 12.5 deg 2 map, Planck 
should be able to detect with a high significance a non-Gaussian level as weak as 10% 
in temperature standard deviation (rms) (5% in Ce), whereas a marginal detection 
would be possible for MAP with a non-Gaussian level around 30% in temperature 
(15% in C e ). 



1 INTRODUCTION 



From the now well measured temperature fluctuations in 
the cosmological microwave background (CMB) we can gain 
some unvaluable constraints on the physics of the early 
universe. An important issue lies in determining whether 
these fluctuations are of Gaussian nature or not. Indeed, 
most of inflation scenarios predict a very low level of non- 
Gaussianity while models involving topological defects can 
give rise to significantly non Gaussian fluctuations. Even if 
a high degree of non-Gaussianity seems disfavored by re- 



For example, for a sufficiently smooth and non degenerate 
random field it is possible to measure the Euler characteris- 
tic (or genus) of the excursion sets via a counting of critical 
points^, classified according to their curvature, i.e. max- 
ima, minima or saddle points in our 2D case (e.g. Adler 
1981, p. 87). Critical points counting is a well known prac- 
tical issue in the context of CMB analysis and specific pre- 
dictions can be derived for smooth random Gaussian fields 
(e.g., Adler 1981), which allow one to constrain the degree of 
non Gaussiani ty of CMB maps by measuring e.g . the Euler 



cents measurements ( Ncttcrficld et al. 2 001 ;_ Lee et al. 2001 



Halverson et al. 2001; Santos et al. 2001: Polenta et al. 2002 



Gott et al. 199C 



Luo 1994 



characteristic (Coles et al. 1989 ; 
Smoot et al. 1994; iBarrciro et al. 200l[) or pe a k statistics 
(Bond and Efstathiou 1987); |Coles et al 1989] |Novikov & 



), the presence of topological defects at a significant level is 
not yet ruled out by observations, as advocated recently by 
e.g. Bouchet et al. (2000), who considered a mixed model 
where Cosmic Microwave Background (CMB) temperature 
fluctuations are seeded in part by cosmic strings. 

There are numerous ways of probing non Gaussian fea- 
tures of a random field such as a CMB temperature fluctu- 
ations map. Since statistical properties of random Gaussian 
fields are entirely determined by their two-point correlation 
function, a natural approach consists in measuring higher 
order correlation functions or related statistics such as cu- 
mulants of the distribution. The measurements can be done 



Jorgcnscn 1996 ; Heavens fe Sheth 1999| ; |Heavens fe Gupta 



1995 



Ferreira et 



al. 1998; 


Magueijo 2000; |Banday et al. 2000|; Gangui & Mar- 


tin 2000a 


; Gangui & Martin 2000b; Verde & Heavens 200ll 


Komatsu et al. 2001; Santos et al. 2001; Kunz et al. 200l|) or 



Forni 1999 



in the space of wavelet coefficients (Aghanim 
Aghanim fe Forni 200l| ). 

Alternatively, non Gaussian features of CMB maps can 
be probed by the topological analysis their excursions sets, 
defined for a random field T(6) as A U (T) = {9 j T{9) > it}. 



200 1| ) 

Following an approach advocated recently in the large 
scale structure context (3D) by Colombi, Pogosyan & 
Souradeep (2000), we propose in this paper to extend the 
counting to ordinary points, i.e. to classify all the points 
according to their local curvature as belonging to "hills", 
"lakes" or "saddles" (see Fig. |l| for an example). By mea- 
suring the relative abundances, Vhm, Piake and Twaddle, of 
these three types of points for various excursions sets and 
smoothing scales, i.e. by exploring the correlation between 
height and curvature^, we are able to extract a mathemat- 
ically well defined Gaussian signature, depending formally 
on only one specific ratio of spectral parameters. As a re- 
sult, a comparison of the measured abundances with those 
predicted in the Gaussian case allows us to detect a certain 
level of non-Gaussianity. 



* The points where the gradient of the field cancels. 

t In this sense, our work is similar in spirit as in Takada (2001). 
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Figure 1. Example of location of "lake", "hill" and "saddle" points in a Gaussian random field smoothed with a Gaussian window of 
size 5 pixels. The points are colored according to their type: "hill" points are in white, "lake" points in dark grey and "saddle" points 
in light grey. Lakes are connected to light grey hills by dark grey saddles. "Lake" and "hill" points are of comparable abundance while 
"saddle" points are the most common. The evolution of the relative abundances with thresholding can be well visualized from this plot: 
for example, when the threshold gets higher, the relative abundance of lake points decreases, while on the contrary the fraction of hill 
points increases, as expected. In this paper we shall rely on details of these variations to discriminate between a Gaussian and a non 
Gaussian random field. 



This paper is organized as follows. In section § g, af- 
ter recalling some useful results about 2D Gaussian random 
fields, we derive analytic predictions for the abundances, by 
extending the work of Bond & Efstathiou (1987). In sec- 
tion § ^| we discuss technical issues involved in the practical 
measurements of 7\m, T'lake an d 'Psaddie- In section § ^ we 
build a x 2 statistic that allows us to combine a set of mea- 
surements in order to quantify the likelihood of an image to 
be Gaussian, from the standpoint of this measure. We then 
apply the method to the case of noisy mixed models involv- 
ing cosmic strings plus cold dark matter. Finally, results are 
discussed in section 8 H. 



2 THEORETICAL CONSIDERATIONS 

In this section we recall some specific properties of 2D 
Gaussian random fields, following Bond & Efstathiou (1987) 
(BE87) (a 2D transcription of the 3D formalism of Bardeen 
et al. 1986). We extend the calculations of this work to ob- 
tain theoretical expressions for the abundances of interest, 
i.e. the fractions of space occupied respectively by hill, lake 
and saddle regions. 



2.1 Facts about a Gaussian random field 

Let us first consider a temperature fluctuation field, 8t(y) = 
T(r)/T — 1 whose 2D Fourier transform is 5t(1<) = 
exp(ik.r) and whose power spectrum P(k) is 



(<5 T (k)<5 T (k')) = (27r) 2 P(k)fo(k + k') 



where 5rj(k) is the Dirac distribution. We can define the 
moments of the power spectrum 



and the ratio 



tf k P(k)k 2 



7 : 



OW2 



(2) 



(3) 



Note that aa is the rms of the fluctuation field St. Natu- 
rally, in the flat sky approximation, i.e. large multipole I 
and small separation angles, the Fourier (flat) power spec- 
trum is related to spherical harmonic power spectrum as 
follows 



P(k) ~ C t and k : 



(4) 



Let us consider from now on only the normalized fluc- 
tuation field 



5(x) = 



(Jo 



(5) 



At a given point r, we note the gradient of the field V5 = r\ 
and the Hessian matrix, = dS/dxidxj. This symmetric 
real matrix can be diagonalised by applying a rotation of an 
angle 9 and we note Ai and A2 the opposite of its eigenvalues, 
with Ai > A2. The normalized trace of the curvature matrix 
reads, 



(1) 



Ai + A 2 
x = , 

and the ellipticity e is defined as 



(6) 
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A1-A2 () 
2a 2 x V ; 

Note that with this notation, x and e should have the same 
sign. If A 2 > then x > and < e < 1/2, if Ai < then 
x < and —1/2 < e < 0. If Ai and A2 have opposite signs, 
then neither x nor e are restricted. 

For a Gaussian random field, the joint probability dis- 
tribution function of 5, rj, x, e and 9 is given by [Eq. (A1.6) 
in BE87] 

V(5,rj,x,e,6)dS dx de dd d 2 r\ = 

s 2 /2 ^d8 ^ {x ^ 5 f/2 dx 
e — e — x 



-r, 2 /a 2 d?rj_ -4a 



2 .2 



8x 2 ede — , 

7T 



where we set 
M=(l-7 2 ) 



'1/2 



(8) 



(9) 



Knowing this probability distribution function, it is easy to 
show that (S 2 ) — (x 2 ) — {rj 2 ) — 1, that 77 is correlated with 
neither 5 nor x nor e, but that (dx) = 7, i.e. the height of a 
point is correlated with the curvature at this point. This lat- 
test result is of particular interest to us since the so-induced 
non-trivial dependence on thresholding will generate a spe- 
cific gaussian signature that we will exploit in the following. 

2.2 Studying the local curvature 

BE87 calculated the density probability of extrema for a 
Gaussian random field as a function of threshold. So using 
the properti es of first derivatives, they open the way to fur - 



ther works (Heavens fc Sheth 1999; Heavens fc Gupta 2001) 



that illustrate the use of the number and spatial distribution 
of extrema to characterize Gaussian random fields. 

We aim at extending this work to second derivatives by 
characterizing the local curvature at any point and then by 
comparing the abundance of 3 defined types. 

The local curvature is defined by the Hessian. By con- 
sidering the sign of its eigenvalues (— Ai and — A2) we are 
led to distinguish three families of points : 

• the "hill" points for which both eigenvalues are nega- 
tive, i.e. x > and < e < 1/2; 

• the "lake" points for which both eigenvalues are posi- 
tive, i.e. x < and —1/2 < e < 0; 

• the "saddle" points for which the eigenvalues are of dif- 
ferent signs. 

The first, second and third family incorporate respectively 
maxima, minima and saddle points. 

Extending the calculation of BE87 we now compute the 
probability that a point (5, r/, x, e, 9) above a given thresh- 
old 5th, i.e. such that 5 > 5th, belongs to any of these 
three classes. More precisely we calculate the quantities 
Pwu(<5th) ee V(0 < x,0 < e < l/2|5 th < S), Pi ake (5th) ee 
T(0 > x,0 > e > l/2|5 th < 5) and P sad die(5th) = 1 - 

■Phill(5th) — ^Plakc(5th)- 

To compute the quantity 'Phiii(5th) we first estimate 
the probability P(0 < x, < e < 1/2, 6 th < S). Starting 
from the distribution, Eq. (^), for the cases of interest to 
us, we can perform straightforwardly both the integration 



I^lo d 2 T]/{' na i) an d J d6/n. Doing so yields a probabil- 
ity function ViS, x,e). Note that the fact that 9 G [0,7r] is 
due to the chosen ordering of the eigenvalues (Ai > A2) im- 
plicit in the distribution Eq. (H) . The subsequent integration 



J de yields the differential density 
A/"hiii(5,K,0 < e < 1/2) dSdx = 



-5 2 /2 



jidS 



-IJ. 2 ( X -^S) 2 /2 



(1 



dx 



(10) 



(2tt) 1 /2 ^ y (27r)V2- 

Then we can still perform analytically the integration over 
x > and we get 



-/Vhm(5, < x, < e < l/2)dS = 
-s 2 /2 ydS 



2V2-K 

- M V<5 2 /(2+M 2 ) 



ill- 



(^1+Erf( 



2J2- 



■ M 



(11) 



The integration over some threshold 5 t h cannot be per- 
formed analytically. However, knowing that 



P(«Sth < S) = -Erfc 



5th 
V2 



it is easy to evaluate numerically the quantity 
f s °° A/"hiii(5, < x, < e < l/2)d5 



Phm(5th) = 



V(6 th < 5) 



(12) 



(13) 



We can still obtain analytically the limiting value "Phm(5th = 
—00), i.e. the fraction of "hill" points in the absence of 
threshold. We find 



7\m(5th = -00) = - M - — 



0.2113 



(14) 



An analogous calculation can be performed for "lake" 
points. It leads to the differential density 



AW(5,:r < 0, -1/2 < e < 0)d5 
-s 2 /2 Hd8 



2V2tv 

- M V<5 2 /(2+M 2 ) 



/' 



-Erfc 



V2 
7M 2 5 



( 7m 2 \ 



(15) 



from which we can deduce numerically 7 3 i a kc(5th) as in 
Eq. (|l3|). The asymptotic value - Pi a i cc (5th = —00) can also be 
obtained analytically, and from parity argument one finds 



■Pl a ke(5th = —OO) 



Phiii(5 t 



■A 



■ — 00) 

) ~ 0.2113. 



(16) 



Two remarkable properties of Gaussian random fields 
emerge from the above analytical calculations: 

• the evolution of the "hill", "lake" and "saddle" point 
fractions as functions of threshold 5 t h depends only on the 
spectral parameter 7; 

• the asymptotic value for 5 t h — » —00 is independent of 
any spectral parameter, i.e. is the same for any Gaussian 
random field. This would not be true if we were considering 
only maxima, minima or saddle points (see BE87). 

To illustrate these results, we draw on Fig. |i| the func- 
tions "Pi a ke(5th) and "Phm^th) for a set of 7 values. Except in 
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Figure 2. Functions 'Phill(^th) an d "Piake(<5th) in the Gaussian 
case for various values of 7 as indicated on the figure. Their value 
is independent of 7 in the limit <5 th — * —00, but their evolution 
with is rather sensitive to 7. Similar conclusions hold for func- 
tion P sad dic(5th) = l~'Phiii(<5th)-'Piakc(<5th), which is not plotted 
for simplicity. 



the limit 8 t h — ► —00, the evolution of these functions with 
5th is rather sensitive to 7. 

At this point, it is important to be aware that till now, 
we have supposed an idealistic, infinite resolution experi- 
ment, i.e. we did not take into account any beam smearing 
effect that would affect any real measurement. Consequently, 
in order to compare the predictions to any true measure- 
ments, it would be necessary to incorporate this effect in 
our calculations. To do so, we should consider the convolved 
field, S * B, where B stands for the instrumental beam re- 
sponse, instead of the idealistic fluctuation field, S. Since 
linear transformations such as the convolution by an arbi- 
trary beam function do not change the Gaussian nature of 
a random field, the beam smearing should appear only as a 
dependence of 7 on the smoothing scale <Tb, that we will note 
7(<Tb). Furthermore, although the problem of beam convo- 
lution in its full generality is very intricate, it turns out to 
be analytically tractable for a Gaussian, symmetric beam, a 
reasonable approximation in practice. We will illustrate this 
point in more details below. 



Figure 3. Measurement of functions Phill an d ^lake 011 Gaus- 
sian simulated maps in the standard CDM case. The trian- 
gles and the losanges give respectively the hill and lake point 
fractions obtained from an average over 200 realizations of size 
12.5 X 12.5 deg 2 . The error bars on each symbol are obtained 
from the dispersion over the 200 realizations. Prior to measure- 
ment, a smoothing with a Gaussian window of size o^, = 50 p i x 
was performed. The continuous and dashed curves correspond 
to the theoretical expectations for a Gaussian random field with 
7cff = 0.48. 



simulated maps. We henceforth limit ourselves to small sim- 
ulated square patches of the sky of width ~ 12.5 deg. These 
patches are considered as being flat and pixelized with a 
512 x 512 Cartesian mesh of resolution # p j x = 1.5'. 

This section is organized as follows: we first describe 
the measurement principles (§ |j.l| ) and apply them to noise 
free realization of Cold Dark Matter (CDM) (Q m = 0.3, 
JIa = 0.7, h = 0.5, Qt = 0.05 and n = 1) Gaussian maps 
generated from their power spectrum (taking into account 
both the uniform distribution of the phases of the Fourier 
m odes and the Rayleigh distribution of their modulus) ^] 
(§ 3.2); in particular, we examine the Gaussianity of the 
distribution function of measurements in order to be able 
to perform later \ 2 tests; finally we discuss practical i ssue s 
related to finite resolution and finite volume effects (§ 3.3). 



3 CONFRONTING PREDICTIONS AND 
MEASUREMENTS ON SIMULATIONS 

In order to test our ability to measure the functions 7\ili(<5th) 
and Piakc(5th)J3 and to measure the incertainties for maps 
of limited extent, we performed several measurements on 



3.1 Principles of measurements 

The measurements consist in determining the fractions 7\ai 
and "Piake for an ensemble of excursion sets (subsets for which 
5 > <5th) of an image smoothed at different scales <7b- To 
insure sufficient differentiability, i.e. at least continuity of 



t Again, we restrict here our analysis to these two functions, since 
function Twaddle (<5th)i which is equal (by construction) to 1 — 
■Phill(<5th) — T'lakeC^th)! does not yield any further information. 



° The power spectrum has been compute d making use 
of the public ly available C MBFAST code (jSeliak fc Zaf 



darriaga 199fi| ) available at http://physics.nyu.edu/matiasz/ 
CMBFAST/cmbfast.html 
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second order derivatives, we choose a Gaussian smoothing 
window. 

Similarly as in Colombi et al. (2000) for the 3D case, 
to measure the local curvature at a given point, we find 
the quadratic form which fits this point and its 8 closest 
neighbors. We deduce from the quadratic form coefficients 
the local values of the second derivatives, i.e. the Hessian 
matrix, then diagonalize this matrix and determine the sign 
of its eigenvalues Ai and A2. Both easy to implement and 
fast, this method turns out to be quite robust in the presence 
of noise. 



3.2 Example: application to noise free realizations 

As a first example, we apply this measurement technique 
to 200 simulated noise free realizations of CDM like tem- 
perature maps smoothed by a Gaussian beam of <7b = 
50 p ix = 7.5'. In each realization, we measure both 7\ili(^th) 
and "Fiake (<5th) for the ensemble of excursion sets defined by 
5th = -3.0, -2.5, -2.0 . . . 3.0, as illustrated by Fig. [| For 
comparison purpose, we plotted on this figure the theoretical 
curves corresponding to the expected value of the parameter 
7(<7b = 5') = 0.48. This effective 7 can be easily calcu- 
lated for a Gaussian beam, since in this context we just 
have to replace P(k) (an exact theoretical input here) by 
P (k) exp(— fc 2 Qh) in Eqs. ^ and (^), (see Fig. g and section 



3.3 for a mode detailled illustration of this calculation). 



With this choice of 7, the theoretical curves fit very well the 
data. The dispersion over the measurements is very tight ex- 
cept at high threshold where we enter a rare events regime, 
as indicated by the larger la error bars. 

To illustrate more visually the protocol, Fig. [| displays 
one CDM realization (noise free again) smoothed at <7b = 



30 P 



the corresponding excursion set for 5 > 1 as well as 



the local curvature map, where hill, lake and saddle points 
are identified by respectively black, white and grey points. 
A cross examination of these images provides some insight 
on the processe involved, i. e. the correlation of extrema and 
local curvature. 

To further investigate the dispersion over the measured 
values of "Phm(5th) at a given 5th, we drawn on Fig. |]] 
the probability distribution function of the hill point frac- 
tion measured from the 200 realizations [similar results can 
be obtained for 7 3 i a ke(5th)]- Two thresholds are considered, 
5 t h = —2 (left pannel) and 5 t h = (right pannel). The 
dashed line on each pannel corresponds to the Gaussian limit 
with same standard deviation as the one actually measured. 
It fits qualitatively well the measurements, so the previously 
drawn 1 a errors are meaningful estimates of the cosmic vari- 
ance errors]^] We checked wether the distribution function 
of the measurements is Gaussian for other values of 5th and 
found that it is indeed the case except at high threshold, 
5 t h > 2.5, where one enters the rare events regime. 

Thus, the measurements in this idealistic noise free case 
are in very good agreement with theoretical expectations. 
Furthermore, the cosmic errors associated with these mea- 
surements are very tight for a Gaussian field. This last point 

^ However, one must keep in mind that error bars for different 
values of 8 t ^ are correlated with each other, but this will be taken 
into account in the analyses conducted later. 




Figure 4. Distribution function of hill points fraction obtained 
from 200 CDM noise free realizations of temperature maps 
smoothed with a gaussian window of size <j b = 30 p i x , similarly as 
in Fig. ^. Two excursion sets are considered, one with <5 t h = —2 
(left panel) and the other one with <5 t h = (right panel). The 
smooth dotted-dashed curve on each panel corresponds to a Gaus- 
sian of same variance as the distribution function. It superposes 
well to the measurements (histograms) for the excursion sets con- 
sidered here. 



will be of particular interest to us since it makes the evolu- 
tion of functions "Phm(5th) and 7 ? i a k c (5th) with 5 t h a sharp 
signature of random Gaussian fields. 



3.3 Measuring 7: spurious effects and available 
dynamic range 

In the above discussion, we used a predicted value of 7 to 
check if the measured "hill" , "saddle" and "lake" point frac- 
tions agreed with theoretical predictions for a Gaussian ran- 
dom field. Therefore, the only thing we proved so far is that 
for a finite realization of a Gaussian random field smoothed 
with a Gaussian of width <7b = 36* p i x , whose (unsmoothed) 
power spectrum is perfectly known, the predicted value 7 
agrees very well with the measured points. In practice how- 
ever, the knowledge of the power spectrum might be less ac- 
curate and we want to develop a non-Gaussianity test which 
can be performed both independently or in combination 
with the power spectrum measurement. Thus it is impor- 
tant to determine from the measurements of these fractions 
themselves an effective value, "/ c g, that fits well the data. 
We will see that for a finite realization of a Gaussian ran- 
dom field, there always exists a value 7 e ft which is such that 
measurements agree very well with analytic predictions. 

As we shall see below, the determination of the 7 e fl 
value is easy and this fact, by itself, is highly significant and 
is the key point of our paper: for a map in which temper- 
ature fluctuations are partly seeded by cosmic strings, we 
shall see that it is actually not possible to find a value of 
7eH leading to point fractions matching the measurements 
in all the available dynamic range - namely, for all possible 
values of 5 t h and at,- However, as we shall see later, it is 
not needed to accurately determine the real value of 7 to 
efficiently constrain the level of non Gaussianity of a map. 

However, studying the difference between 7 and r y e s can 
help to determine the available dynamic range, i.e. the set 
of values of 5th and ah in which one can trust the mea- 
surements. We already noticed in previous section that the 
density threshold should be small enough, 5 t h ^ 2.5, to avoid 
entering the rare event regime, where the cosmic distribu- 
tion function of measurements becomes non Gaussian. Here, 
we are concerned by two effects that we ignored previously: 

• Pixelization effects, or equivalently, effects of finite reso- 
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Figure 5. Illustration of the measurement protocol. Considering a CDM temperature map with 512 X 512 pixels, of width 12.5 deg X 
12.5 deg and smoothed with a Gaussian window of size 3 pixels (top left panel), we identify "hill", "lake" and "saddle" points according 
to their local curvature (respectively by dark grey, white and light grey points on the bottom left panel). The fraction of space occupied 
by each kind of point is studied as a function of the density threshold. As an example, the excursion set corresponding to 8^ = 1 is 
shown in right pannels, which are the same than left panels except that only regions with <5 > <5 t h are shown, the rest being coded in 
black. As one can see, the relative abundance of "lake" points in bottom right panel is now smaller than on bottom left panel, while 
which of "hill" points has augmented, as expected. The measurement of the relative fraction of the three kind of points as a function of 
density threshold is the key idea of this paper. 



lution: as discussed above, it is important to smooth the data 
to insure sufficient differentiability. The smoothing scale 
should be large enough compared to the pixel size to avoid 
anisotropies or discreteness effects brought by the pixeliza- 
tion. 

• Finite volume and edge effects: our experimental maps 
cover a rather small fraction of the sky, due to the lim- 
i ted dynamical range in the cosmic-string simu lations from 



lnPhiii(5th 



b-< 



(17) 



(Bouchet et al. 198c; Bennett & Bouchet 199C). Therefore 



to avoid reducing too much the number of statistically inde- 
pendent regions of the map, the smoothing scale should not 
be too large. In a more realistic experiment covering a large 
fraction of the sky, finite volume and edge effects should not 
be as much as of a concern. 

The practical measurement of y e g is made straight- 
foward by noticing that the dependence of the function 
Vhi\\(Sth = 1) on 7 is very well approximated by an exponen- 
tial law|^in the domain of interest to us, i.e. 0.4 < 7 < 0.95: 
the fit 



II Note that such a fit is also appropriate for values of <5 t h different 
from 1. 



with a = -1.4099 and 6 = 1.2395 is accurate to 0.6%. The 
choice of the particular value <5 t h = 1 is an ad hoc compro- 
mise coming from the competition between two effects: (i) 
the dependence on 7 of the quantity Pull (5th) increases with 
5th (Fig- 12|), but (ii) the uncertainty on the determination of 
Vhi\\(Sth) increases with 8th (Fig. g). 

To test the effects of finite coverage and finite resolution, 
we examine a CDM case, where we generated once again 
some Gaussian field realisations from their Ce- Assuming 
as before that the field has been smoothed by a Gaussian 
of width <Tb, we can easily compute the prediction for the 
function 7(<Xb). 

In Fig. kj. the measured values of 7 c ff is displayed in 
each case as a function of smoothing scale. An average over 
50 realizations of 512 x 512 pixels maps is performed. The 
error bars on the figure represent the corresponding scatter, 
which increases with <7b as expected, due to finite volume 
effects. The solid line corresponds to the theoretical func- 
tion 7(<Tb)- Agreement between 7 c h and the analytic predic- 
tion is good, except at small scales where pixelisation effects 
contaminate this measurements. We however see that pix- 
elisation effects become negligible when <7b/# P ix > 3. Note 
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0.4 







5 10 15 20 
Smoothing a h / 6 pix 



Figure 6. Measurement of -y e g as a function of the smoothing 
scale Ob/(?pj x . The solid line gives the theoretical values of 7 for 
various o-j,. They were computed using the power spectrum Ci. 
The crosses with error bars correspond to measurements of 7 e ff 
using the measurements of 'Phill('5th = 1.0). The error bars cor- 
repond to ±1 the measured rms on 50 realisations. Pixelisation 
effects bias the measurements of -y e g till (7b/0 p i x = 3. The error 
bars grows with as expected from sampling variance argu- 
ments. 



that there should be an upper bound as well for Ob as dis- 
cussed above, since cosmic errors increase with scale. Our \ 2 
analysis below will anyway naturally take that into account 
by giving a lesser weight to scales with larger errors. This 
reasonable agreement between the measured 7<,g at a given 
scale Ob and the expected 7(ob) validates the approach we 
will use in the next section which consists in finding from 
the fraction themselves an ad hoc j e g, and deducing from it 
an agreement with Gaussianity. 



4 TESTING NON-GAUSSIANITY : 

PRINCIPLE AND APPLICATION TO 
MIXED MODELS 

^,From the previous section, we conclude that in the case of 
a smoothed Gaussian random field, the functions Phui(<5th) 
and 'Pvoid^th) can be measured accurately and fit very well 
the analytical predictions provided 7 be considered as an ad- 
justable parameter. Furthermore, their probability distribu- 
tion function is well approximated by a Gaussian if 5 t h ^> 2.5, 
which now allows us to define a rigorous measurement pro- 
tocol based on \ analysis that we shall apply to simulated 
data in the trustable scale range determined above, namely 
o- b > 3. 

In this section, we first detail how we use the functions 
7\iii(<5th) and "Pvoid(5th) to build a \ 2 statistic testing Gaus- 
sianity (§ [4.l[ ). Then we apply the method to simple simu- 
lated mixed models involving cosmic strings plus CDM and 
see how it can be used to bound qua ntitatively the relative 
contribution of cosmic strings (§ 4.2). 



4.1 The method 

Let us assume as before that we have at our disposal an ob- 
served temperature map St, on which we measured the frac- 



tions "Phiii(5thi, ObJ and "Flake (5th;, ObJ for a set of threshold 
values {^thj} and smoothing scales {obj}- If the power spec- 
trum of this random field is known, we can in principle de- 
fine analytically a set of values of 7 for each smoothing scale, 
{7 CTb }. Assuming furthermore that the field is Gaussian and 
that our estimators of "P h iii and "Pi a k c are unbiased, it is easy 
to deduce from this 7 set some theoretical expectations for 
the set, VxiSthi , ObJ = ("Px^ttj , 0"bi)) , where X stands from 
now on for "lake" or "hill" . To quantify the distance between 
these theoretical predictions and the measurements, we in- 
troduce the standard \ 2 statistic. Since these measurements 
can obviously not be considered as independent, we have to 
introduce the full theoretical variance-covariance matrix 



Cr 



({Vi-Vi){Vv-Vv)), 



(18) 



where we use the short hand notation Vi = "Px(5th ; , Cbj) and 
Vy — "Px' (^thj, , ob., ) • Then, the statistic 



E 

IF 



(V1-V1) (C n {Vv -Vv) 



(19) 



is expected to follow a x -distribution, as illu strated by a 
practical example below. Indeed results of § 3.2 suggest that 



the distribution function of the measured V\ is nearly Gaus- 
sian if 5th ^ 2.5. 

In practice, we have to introduce two more subtleties in 
this protocol. 

First, as discussed extensively in 



3.3 



the practical re- 
alization of a Gaussian random field yields a measured func- 
tion V\ matching the theory, V, but with an effective value 
of 7. Taking into account some additional Gaussian noise, 
as discussed in more details below, would make this effect 
even stronger. Rigorously, we should modelize this 7 c ff effect 
in terms of statistical bias on the estimators of the function 
Px^thj °~b), but we proceed here differently, for simplicity: 
we extract the 7 e ff(o"b) from the data by fitting the analytic 
prediction for function "Phm(<5th = l,°"b) to the measured 
one. Thus, our "theory" depends itself on the data through 
7off. As already explained, the critical test for non Gaus- 
sianity will be on the evolution of the functions Vx{o"th, o~b) 
with <5 t h, which, once 7 c ff is determined, is entirely fixed. The 
practical measurement of 7 e ff is performed as explained in 
the previous section. In principle we could include the value 
of 5th used to measure 7 c ff as a varying parameter in our \ 2 
test, to optimize the analysis. Given the level of accuracy 
of our numerical experiments and since we just aim here to 
illustrate the method in a simple and convincing way, we 
did not feel necessary to do so. 

Second, even if in principle we could compute ana- 
lytically or numerically the variance-covariance matrix, for 
the sake of simplicity we evaluated it using a Monte-Carlo 
method based on a few hundreds of Gaussian realizations 
having the same power-spectrum and noise properties as the 
input map. So for each model we will consider below (pure 
CDM, mixed model with CDM + cosmic strings, and pure 
cosmic strings), the covariance matrix C should be differ- 
ent in each case. However, in practice we checked that this 
matrix coefficients depend very slightly on 7 e g in the range 
of interest to us, i.e. 0.4 < 7 < 0.95 so that we consider 
for simplicity only the same C matrix in each case. Another 
issue concerning this matrix is that it might be singular. 
Indeed, since fractions measured at various thresholds and 
different scales can be highly correlated, e.g. the values at 
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very low threshold where the dependence on 7 is very weak, 
the studied {5th;} and {cn^} sets have to be restricted so 
that C is invertible. 



4.2 Application on mixed models for MAP and 
Planck 

As an illustration of the accuracy of our statistics, we apply 
it in the framework of MAP^ and Planck Surveyor^ ex- 
periments to a mixed model where temperature fluctations 
are seeded in part by cosmic strings and otherwise by adi- 
abatic inflationary perturbations in standard CDM model 
{e.g. Bouchet et al. 2001 and references therein). Our ap- 
proach is rather simple, since we add together a CDM map 
and a temperature fluctuation map obained from ray-tracing 



in cosmic string simulations (Bouchet et al. 1988 



Bennett 



& Bouchet 199C) and neglect possible cross-correlations be- 
tween such maps. This implicitely supposes the existence 
of a standard CDM inflationary epoch followed by a phase 
transition during which cosmic strings appear and then im- 
print supplementary fluctuations on the dark matter distri- 
bution. Neglecting cross-correlations mentionned just above 
is equivalent to assume that the distribution of cosmic string 
as well as their dynamical evolution is not coupled with that 
of cold dark matter, which should be a good assumpt ion at 
the level of approximation considere d in this paper (Linde 
fc Riotto 19971 |Contaldi et al. 200l| ). We neglect also any 
contribution to the fluctuations prior to the last scatter- 
ing surface (lss) so that the comparaison with results from 
Bouchet et al. 2001, who considered this contribution, is not 
immediate at small scales. 

To take into account the noise expected in the best 
channels of MAP and Planck Surveyor experiments (respec- 
tively 12.8 /iK/K per 0.3 x 0.3 deg 2 pixels and 2 fiK/K per 
8' x 8' deg 2 ), we add an extra Gaussian white noise n to 
our square maps. Note again, that as long as the noise is 
Gaussian, its presence can be simply taken into account by 
a change in the value of 7 c ff. Thus, our simulated tempera- 
ture maps read 



(1 - 0) S 



strin; 
T 



. a r CDM 

+ p S T + n, 



(20) 



where 1 — represents the fraction of the signal seeded by 



cosmic strings, given the fact that we impose rms(<5, 



:tring 



) = 



rms((5£ DM ) = rms(^ OBE ) for the fields smoothed at a 
scale corresponding to the COBE-DMR beam, i.e. <Tb = 
7°/V81n2. Note that we differ from Bouchet et al. 2001 
in the relative normalisation of the two contributions since 
they define their relative normalisation with Ce, i.e. they 
considered a family of models parametrized by Ce = (1 



„ \ .^string . y^fL 

a) C e + a C? 



where both C\ 



string 



and Cf DM are 



COBE normalised and where C| trmg includes contributions 
from before and after the lss. If we wanted however to ex- 
press our results, in terms of 0, we should consider 1 — a = 
(l-/3) 2 /((l-/3) 2 +/3 2 ). 

We now construct the previously defined x 2 statistics 
for various values of and try to figure out what is the 
smallest "string like contribution" that our method should 
be able to detect with a high statistical significance. 

http://map.gsfc.nasa.gov 
tt http : //astro . estec . esa . nl/SA-general/Proj ects/Planck/ 
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Figure 7. Histogram of the measured probability distribution 
of the quantity \ 2 defined by Eq. (|is|), from 300 noisy (Planck 
best channel noise level) realizations of CDM temperature fluc- 
tuations. The dot-dashed curve corresponds to a x 2 -distribution 
with 47 degrees of freedom (dashed line). The good agreement 
allows us to use x 2 as a measure of the Gaussian nature of the 
observed field. 



4-2.1 Constructing the x 2 and adequacy of the \ 2 
probability distribution 

Even if one would like to include as many threshold values 
and smoothing scales as possible in the analysis, this is in 
practice not necessary. Indeed, some measurements are, at 
least in the Gaussian case, highly correlated with each other, 
due to e.g. some strong constraints (in the very low threshold 
regime, 8 t h —* —00, both 7 3 iake(<5th) and 'Phm(<5th) tend to an 
asymptotic value independent of the spectral parameters). 
Therefore, one can restrict the set of values of (o"B,5 t h) by 
requiring that the C matrix (see Eq.[l8]) be non-singular. 

By investigating the numerical properties of C, we 
find for all the maps q| that <5 t h = —2,-1.5 . . ., +1.5 and 
0"b/# P ix = 6, 10 and 14 (corresponding to repectively 21', 
35' and 49' on the sky) are satisfying values for the Planck 
case, whereas 5th — —1.5,-1.0..., +1.5 and <7b/# p ix = 8 
and 12 (corresponding to respectively 28' and 42') are rea- 
sonable for MAP. Of course, this choice of the set of values 
of (o\b, <5th) is strongly influenced by the fact that our map 
is rather small: for a full sky survey, the sampling would 
probably be different. 

Given these restrictions, it is now important to explore 
the properties of the probability distribution function of the 
X 2 statistics defined in Eq. (119). As stated previously, we 
expect it to follow a x 2 -distribution, at least in the Gaus- 
sian case, since the distribution of measured values of Vis 
is well described by a Gaussian (see 



3.1). We tested this 



by generating 300 realizations of noisy (Planck level) CDM 
like only temperature fluctuations for which we measure x* ' ■ 
Its probability distribution is drawn on Fig. [5], where we su- 
perimposed a theoretical x 2 -distribution with 47 degrees of 
freedom (dof) (48 measurements: 6x4 measured hill points 



H As already stated, we checked that this choice is adequate for 
various C matrices, i.e. corresponding to different 7s, altough we 
use only one in practice. 
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Figure 8. Measurement of hill and lake point fractions in a 
Planck simulation with pure CDM (/3 = 1.0). The left column 
shows the measured evolution of 'Pi a i <c as a function of thresh- 
old, 5th, for two smoothing scales: top left panel corresponds to 
°"b/f9pix = 10 and bottom left panel to u h /8 pix = 14. Symbols 
without error bars represent the measurements in the realization, 
while symbols with errorbars correspond to the Gaussian predic- 
tion match ing the value of 'PhinC^th = 1) as explained more in 
detail in § 4.1 The error bars are obtained from the 1-cr disper- 



sion over 300 realizations. The right column of panels is similar 
to left column, but for Phill- 



fraction and 6x4 measured void point fractions; 1 param- 
eter, 7). The agreement is obviously very good, comforting 
the analysis procedure and making this \ 2 statistic adequate 
for assessing Non-Gaussianity.^ 



4-2.2 Ability to distinguish a "Non-Gaussianity Level" 
with a realistic noise level 

Making use of this \ 2 statistics, we are now in position to 
determine how well we can distinguish non-Gaussian signa- 
tures. 

In Figs. ^, ^ and [To], we display some measurements of 
functions 'Piake(o'th) and 'Phiii(o'th) in the Planck case. Three 
values of /3 are considered: pure CDM with /3 — 1, hybrid 
model with (3 = 0.3 and pure string model with /3 — 0. 
Again, On each figure, the Gaussian limit obtained from 
the procedure explained in detail in § ^t] is compared to 
the measurements. The error bars on the theoretical pre- 
dictions correspond to then the lcr dispersion over 300 re- 
alizations of pure random Gaussian fields with same power- 
spectrum and noise properties as the data, as explained in 
end of ! 



4.1 



Very good agreement is found for CDM, as ex- 
pected, while the hydrid model exhibits some discrepancies, 
especially for "Piakc- These discrepancies are overwhelming in 
the string only model. Note that since we fit the measured 



§3 Here, we do not test if the x 2 estimator always follows a \ 2 ~ 
distribution, i.e. also for non-Gaussian maps. Such a test is not 
necessary since the prior we use to analyse the data assumes un- 
derlying Gaussianity (it would be impossible to do otherwise in 
practice). 



Figure 9. Same as Fig. g but for the an hybrid model with 
P = 0.3. 
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Figure 10. Same as Fig. ^ but for the pure string model (fl 
0.0). 



value of "Phm(<5th — 1) to estimate the Gaussian prediction, 
one expects non Gaussian features to show up more in the 
measured shape of function Viake(Sth) than that of function 
7\iii(£'tii), as indeed seen at least for the mixed model. Note 
also that if for the mixed model, the measured "Pi a ko(5th) 
tends to lie below the Gaussian prediction, the reverse hap- 
pens for the pure string model (this seemingly counterintu- 
itive result stems from our comparing the data to a "the- 
oretical" curve which is different in the 2 cases) . Finally, 
it is important as well to recall that data points plotted on 
Figs. P |i"o| represent only 2/3 of all the points used for the 
X 2 computation since we consider one more smoothing scale, 
Ob/^pix = 6. 

To illustrate more qualitatively these measurements, 
Figures |ll] and |l3| display, similarly as in Fig. ^, an example 
of initial smoothed noise free map, its thresholded counter- 
part (with 5 t h = 1) and the corresponding local curvature 
maps with lake, saddle and hill points. Figure [ll] corresponds 
to the hybrid model case (P = 0.3) while Figure |l| corre- 
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Figure 11. Same as in Fig. |j but for a noise free mixed model with /3 = 0.7. Note how difficult it is to distinguish it from the mere 
CDM map. 



sponds to the string only model (/3 = 0.0). A strong simi- 
larity can be seen between the mere CDM and the hybrid 
models despite the differences in measured "Pi a ke functions. 
Note also the peculiar patterns in the pure string curvature 
map: there seem to be rather extended saddle point regions, 
which calls for other specific pattern detection statistics. 

Table |E| summarizes the results obtained for x 2 using 
the noise level of the best channel of either Planck or MAP. 
It shows, for various values of f3, the measured x 2 > the \ 2 P er 
dof, as well as the probability (noted Prob.) that a random 
variable, in a ^-distribution with the relevant (number of) 
dof (shown above to be adequate) is greater than the ob- 
tained % 2 - This quantity gives therefore the significance of 
the detection. 

Roughly speaking, the smaller the \ 2 y the better is 
the agreement with a Gaussian distribution. Reversely, the 
higher the obtained \ 2 1 the more unlikely does the analyzed 
map follow a Gaussian distribution. 

Note that the numbers displayed in Table ^ represent 
in no way the best Planck and MAP can provide with our 
method, since we only considered 0.3% of the sky for one 
single channel. Nevertheless, these figures give some trends 
of the ability of our method, the purpose of this paper. 

With these limitations in mind, we see that MAP can 
get a meaningful detection (1% level) of the non-Gaussianity 
induced by string like background with j3 — 0.7, correspond- 
ing to a number 1 — a — 0.15, a number slightly lower than 
the constrains from Bouchet et al. (2000), 1 — a = 0.18. 
Any smaller non-Gaussian contribution is not detected with 



a good significance. On the other hand, we see that Planck 
could obtain, using this method, a highly significant detec- 
tion level, i.e. smaller than 0.1%, for /3 between 0.8 and 0.9, 
i.e. 0.05 < 1 - a < 0.12. 



5 DISCUSSION AND CONCLUSION 

In this paper, we have tested a statistic based on local cur- 
vature measurements by counting three well defined types of 
points according to the sign of the eigenvalues of the Hessian 
of a temperature map. In the first section, we computed the 
relevant theoretical expectations. We validated them and 
demonstrated our ability to properly measure them in the 
next section and eventually introduced this statistic as an 
accurate non-Gaussianity test (by performing a % 2 like anal- 
ysis) in the last sections, where we applied it to the case of 
mixed models with realistic noise levels. Even if our phys- 
ical simulation might be somewhat simplistic (e.g. no pre- 
lss string contribution, no foregrounds residuals, etc.), they 
showed clearly that we are able by using this technique in the 
Planck context, to distinguish a non-Gaussian background 
amplitude around 10% in temperature rms (5% in CV) with 
a good significance using only 0.3% of the sky and one single 
channel. Consequently, these results are definitely positive 
and, as it is, this method seems to be ready to be applied 
on true data. 

Naturally, in practice this test should be used in con- 
junction with other more standard statistical measurements, 
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Figure 13. Same as in Fig. Q but with noise free string induced only temperature fluctuations. Note the interesting structures in the 
curvature map. 







MAP (dof 


= 31) 


Planck (dof 


= 47) 






x 1 


x */dof 


Prob. 


x 2 


X*/dof 


Prob. 


p = 


1.0 


34.1 


1.1 


0.32 


46.6 


0.99 


0.49 


13 = 


0.0 


109.3 


3.52 


3.1 10- 10 


879.3 


18.7 


0.0 


P = 


0.7 


120.9 


3.9 


4.01 


217.1 


4.6 


0.0 


P = 


0.8 


68.2 


2.2 


1.3 10~ 4 


180.5 


3.84 


0.0 


= 


0.9 


55.8 


1.8 


4.0 10~ 3 


101.1 


2.15 


7.0 10~ 6 



Table 1. The measured x 2 > X 2 /dof and significance probability (Prob.) are shown for various values of f) and for our virtual MAP and 
Planck experiments, as explained in the text. 



like 2-, 3-, 4-, ...points correlation functions or their har- 
monic transform, power spectrum, bispectrum, trispectrum 
... or even already presented more sophisticated ones like 
the distribution of wavelet coefficients, etc. This array of 
methods might actually help in identifying the source of 
non-gaussianity detected in a realistic context. Indeed, in 
practical situations, important non-Gaussian effects might 
be induced by CMB contaminants, in particular galactic 
foregrounds, residual map stripping or point sources, or any 
instrumental systematic error. 



this method and so is the pixelisation, provided that the 
smoothing scale is about 3 times the pixel siz ^^| (see §3.3). 

A second advantage of the method lies in the fact that 
it is mathematically well-defined and simple since it depends 
only on one single spectral parameter, 7 that can be either 
measured from the fractions themselves or from the Ct mea- 



surement (see § 3.3). Thus measurements are well understood 
and well controlled from a theoretical point of view. 



5.1 Possible advantages 

One important advantage of these statistic is its local nature, 
i. e. the fact that we are only interested by the 8 neighboring 
pixels. The sphericity of the sky is thus not an issue for 



^ Which is convenient, since pixelisation of low noise experi- 
ments are usually done at a size similar to one-third of the FWHM 
of the beam. 
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5.2 Outlooks 

As we mentioned earlier, by moving to second order deriva- 
tives, this work opens the way to many more characteriza- 
tion of random Gaussian field. We expect that an analogous 
work could be done by considering critical points only, which 
would lead, e.g. , to a measurement of the Euler character- 
istic as a function of threshold. However, we expect this re- 
striction to be less efficient in probing non-Gaussianity since 
it appears that by keeping all types of points we benefit more 
from the full structure of the probability distribution func- 
tion. This point will be demonstrated in a future work. In 
addition, before applying the method to real data, it will also 
be interesting to compare it's ability to that of other esti- 
mators on the same test protocol e, but extended to in clude 
other families of models, as e.g. (Lyth fc Wands 2001). 
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